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Abstract 

The many-body dynamics of a quantum computer can be reduced to 
the time evolution of non-interacting quantum bits in auxiliary fields by 
use of the Hubbard-Stratonovich representation of two-bit quantum gates 
in terms of one-bit gates. This makes it possible to perform the stochastic 
simulation of a quantum algorithm, based on the Monte Carlo evaluation 
of an integral of dimension polynomial in the number of quantum bits. 
As an example, the simulation of the quantum circuit for the Fast Fourier 
Transform is discussed. 



1 Introduction 

The potential use of quantum computers for solving certain classes of problems 
has recently received a considerable amount of attention (see, e.g., jl], || for a 
comprehensive review). Several quantum algorithms have been developed, such 
as quantum factoring Q , having the potential for revolutionizing computer sci- 
ence. The purpose of this paper is to explore the application of a Monte Carlo 
method that has been developed in the context of quantum many-body systems 
to the simulation of quantum computers fri . Quantum computers can be seen 
as peculiar quantum many-body systems that evolve according to a non-local 
time-dependent interaction so as to carry out a "computation" . The component 
quantum bits (qubits) interact via a sequence of quantum gates, performing each 
a prescribed unitary transformation (rotation, Hadamard transformation, con- 
trolled not, controlled phase, etc.) Q. Two-bit (or n-bit) gates therefore effect 
non-local interactions between qubits, and the "quantum algorithm" (charac- 
terized by a network of quantum gates) corresponds to a specific sequence of 
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unitary transformations, i.e., a time- dependent interaction. Numerous methods 
have been developed for years in order to treat general quantum many-particle 
systems (see, e.g., It is therefore intriguing to examine whether the applica- 
tion of the same methods to quantum computers might be similarly successful. 
We focus here on a stochastic approach based on the Hubbard- Stratonovich 
transformation || which has been shown to be suitable for the description of 
quantum many-body systems (see, e.g., @, [l(J). The central idea of this ap- 
proach is to replace the many-body propagator for the entire quantum computer 
(of say L "interacting" qubits) with L one-bit propagators in fluctuating aux- 
iliary fields, thereby "decoupling" the qubits. More specifically, solving the 
quantum dynamics of the L-bit computer in a very high-dimensional Hilbert 
space (d = 2 L ) reduces to evaluating a high-dimensional - but polynomial in 
L - integral over auxiliary fields. The latter is then approximated by use of a 
stochastic method. 

2 Quantum computer as a many-particle system 

Consider a quantum computer consisting of a register of L qubits supplemented 
with a quantum algorithm, defined as a sequence of G quantum gates. The 
total unitary transformation characterizing the quantum computation is thus 
expressed as an ordered product of operators (note the product from right to 
left) 

G 

u = n u g = u G ■ ■ ■ u x (i) 

3=1 

where U g is the unitary transformation performed by the <?-th gate, and G is the 
total number of gates. It has been shown that two-bit gates are universal, that 
is quantum gates operating on one and two qubits are sufficient to construct 
a general quantum circuit ||, [ll], [l^, |T^]. Therefore, we restrict ourselves to 
the simulation of quantum circuits made of G two-bit gates (U g being a two-bit 
gate acting on qubits a g and b g ), keeping in mind that an arbitrary quantum 
computation can be achieved with an appropriate sequence of such gates. (Ob- 
viously, one-bit gates can always be incorporated into two-bit gates.) Note that 
an efficient quantum algorithm must have G polynomial in L. For example, 
the quantum Fast Fourier Transform (FFT) circuit fl4| used in quantum fac- 
toring j| requires G = L(L — l)/2 two-bit gates. A two-bit quantum gate that 
effects a unitary transformation on qubits a g and b g can be written generically 
as the two-bit operator 

U g = e- ia ° A s B <> (2) 

where a g is a real number and A g , B g are two commuting one-bit Hermitian 
operators referring to the qubits involved in the quantum gate [i.e., the operator 
A g {B g ) affects qubit a g (b g )]. For example, the controlled-NOT gate || acting 
on qubit a (as a control) and qubit b (as a target) has a = tt/4, A = 1 — a z , 
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Figure 1: Hubbard-Stratonovich representation of a two-bit quantum gate in 
terms of two one-bit gates in fluctuating auxiliary fields a, t. 
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and B — 1 — a x , with a x , a z being Pauli matrices. The Hubbard-Stratonovich 
representation of U g is obtained by writing the identity 

ICtgAgBg = ia g (Ag - Tg)(Bg ~ <J g) + iOtgOgAg + IttgTgBg ~ Id gO g7 g (3) 

where a g and r g are real auxiliary fields corresponding to the g-th gate, and 
then integrating the exponential of Eq. (|3|) over a g and T g , resulting in 

\ a ol I . , ... ,-• -.- —ia„a„A„ — ia„T„B„ 



Ug = L_9± d(Jg dTg e « 9 < rB T 9e - WsJj A !e -. a „T s i* s (-4) 

This expression is most important because the two-bit gate U g is represented 
as an infinite sum of products of (field-dependent) one-bit gates, e - ia 9 a 9 A 9 an d 
e -ta g T g B g _ p or a gj ven va i ue f th e fields a g , r g , the two qubits a g and b g 
act as non-interacting particles and evolve independently (they do not become 
entangled when initially prepared in a product state). [j] Only the sum over fields 
creates a "coupling" between them, as pictured in Fig. |]. As a consequence, for 
a given set of a g 's and r s 's, the time evolution effected by the whole quantum 
circuit can be computed separately for each qubit: one calculates the time- 
evolution of L qubits in L two-dimensional Hilbert spaces rather than the time- 
evolution of a single quantum state (the state of the entire computer) in the 
full 2 L -dimensional Hilbert space. This exponential reduction of the size of the 
Hilbert space appears clearly when writing the total unitary transformation for 
the quantum circuit 

U = J Da cxp ( iY,<* 9 <r 3 T 9 J \[V g {a g ) Wg(r g ) (5) 



U[a] 



where Da = Y[ g ^rda g dT g is the measure over the auxiliary fields a±, ■ ■ ■ aa and 
r i7 ' ' ' t Gj an d U[a] is the total unitary transformation for a given "path" a in 

1 Several alternative Hubbard-Stratonovich representations of a 2-bit quantum gate requir- 
ing only one auxiliary field per gate can be written, the drawback being that the involved 
1-bit transformations are in general non-unitary. This is under current investigation. 
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auxiliary-field space. Here, V g (a g ) = e~ %CCg<TgA9 and W g (r g ) = e - M s T s s s stand 
for the unitary transformation performed by the one-bit gate acting separately 
on qubit a g and b g , respectively. (The one-bit gates V g and W g replacing the g- 
th two-bit gate U g depend respectively on the auxiliary field a g and r g .) Eq. (^|) 
involves only one-bit operators, and therefore describes the time-evolution of 
L non-interacting qubits (averaged over auxiliary fields). The operator U[o~] is 
more conveniently written as a product of one-bit operators over the L qubits, 

U[a} = f[U^[a} (6) 

1=1 

where the one-bit operator U^ l \ describing the overall evolution of the l-th qubit, 
is expressed as the ordered product 

U^[a] = J{uf{a g ,T g ) (7) 
g 

with 

!V g (a g ) if I = a g 
W g {r g ) if I = b g (8) 
1 otherwise . 

The drawback of this exponential reduction in the Hilbert space is obviously the 
2G-dimensional integral over fields in Eq. (||) , which can only be approximated 
by a numerical method in general. The underlying idea of a stochastic method is 
to compute only the dominant terms in this integral, that is to consider the paths 
in auxiliary- field space that contribute the most to it, assuming that this yields a 
good estimate of the exact integral. Several (more or less efficient) Monte Carlo 
techniques can be thought of for sampling these paths, but a generic "sign" 
problem unavoidably occurs due to the complex weight in Eq. (^). This will be 
discussed later on. However, the central point here is that the dimension 2G of 
this integral is polynomial in the dimension of the problem - i.e., polynomial 
in the number of qubits L - at the condition that G — poly(L). The latter 
condition is fulfilled for any efficient quantum algorithm, suggesting that the 
Monte Carlo simulation of a quantum computer might be interesting if the 
"sign" problem is circumvented. 

3 Stochastic simulation of a quantum computer 

Consider now the stochastic calculation of the quantities of interest in a general 
quantum computation. In the context of quantum many-body systems, stochas- 
tic methods are especially appropriate for calculating quantum expectation val- 
ues, so that our goal is to express the output of the quantum computation as an 
observable. Assume that the entire quantum computer is initially in a product 
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state IO1O2 • • -Ol). (If this is not the case, the first step of the computation 
should simply be the preparation of the correct initial state from IO1O2 • • - Ol).) 
The quantum computation (i.e., the unitary transformation U) is implemented 
by a quantum circuit acting on this initial state. The final step of a quantum 
algorithm is then to measure a set of "output" qubits (not necessarily all the L 
qubits). We restrict ourselves here to quantum algorithms that provide a deter- 
ministic result (unlike quantum factoring). We assume therefore that the output 
bits are in a product state, so they can be measured individually (i.e., one can 
perform an inclusive measurement of each of them separately). The most gen- 
eral observable O with vanishing variance (deterministic output) consists then 
in a product of one-bit observables, and several such O's can be measured si- 
multaneously (If the output bits are not in a product state, one should extend 
the quantum computation with a unitary transformation mapping the entangled 
final state into a product statej^j) More generally, the quantum many-particle 
simulation approach allows us to prescribe the value of certain qubits in the 
output register. We separate the L output qubits into L m measured qubits, L p 
prescribed qubits, and L t — L — L m — L p traced over qubits (i.e., "scratch" 
qubits that are necessary to make the overall computation unitary, but are not 
observed in the final measurement). The observable can then be written as 

O = IJ °m (9) 

{»"} 

where O m is a one-bit observable acting on qubit m. Consequently, the re- 
sult of a quantum computation can be written as the expectation value of the 
observable O, 

/n] = (0i---0l\WOPU\0i---0l) , , 

[U> (0 1 ---0 L \WPU\0 1 ---0 L ) {W) 

where we define the projector P as 

p = U p p ( n ) 

M 

where P p = \tt p ) (ir p \ is a projector on the prescribed value tt p for qubit p. Note 
that it is crucial to consider a quantum algorithm such that the variance of O 
vanishes when the prescribed qubits have the correct value, so that Eq. ( |To| ) 
yields the deterministic output of the quantum computation. The only variance 
in the simulated output will be the statistical noise resulting from the stochastic 
evaluation of Eq. (||) . 

2 It is not clear whether this requirement makes the extended quantum computation much 
harder in a general case. At least, some quantum algorithms are known to provide a deter- 
ministic result, such as Grover's algorithm so that the output bits are then in a product 
state. Note that the same requirement must be met for the recently suggested realization of 
quantum computers using NMR experiments jUl . 
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The central point now is that, using Eqs. (JsJ) and (^), the numerator and 
denominator in Eq. (|l^) can be expressed in terms of an infinite sum of products 
of L one-bit matrix elements for each qubit, 

_ [Do Da' e l ^ s a ^-« ] IltiN^'V] Q [ " P m U [ M \0j> {u) 
f Da D & e'E/.fv.-ir;) UtA°l\ um ^') pll] U[l] \°\ \°l) 

where a' represents the set of auxiliary fields (a' g and r' g ) used in the Hubbard- 
Stratonovich expression of W. This can be written more concisely as 

^ ; fDaDa' exp(-iS[a,a'}) [ ' 

where the (complex) action is defined as 

S[a,a'} = -J2 a a( a 9 T 9 ~ a 'g T 'g) 



+ iJ2MQi\U m [v']P ll] U [l] [<T]\0i) (14) 
1=1 

The operator is a one-bit projector if the i-th qubit is prescribed, and the 
unit operator otherwise. The estimator of O is 

n _{ T (0 l \unt\ [<T >] O m pW[/W[q]|o ; ) 



1=1 



where is the one-bit component of the observable O if the l-th qubit is 
measured, and the unit operator otherwise. Note that the L matrix elements in 
the right-hand side of Eq. ( |l4| ) are for single qubits, so that the calculation of 
the action involves ~ AG products of non-unit 2x2-matrices. (There are 2 fields 
per gate, and the Hermitian conjugate must be considered together with U.) 
The calculation of Eq. ( |i~5| ) requires essentially the same operations. 



4 Sampling of the auxiliary-field paths 

Let us now consider the stochastic evaluation of Eq. (|l3|) based on a sampling 
of the paths (set of cr g 's ,r g 's) that contribute the most to the integral. The 
simplest possibility is to perform an importance sampling of the paths according 
to the weight |e -lS |. (Note that this weight is not equal to 1 since S is generally 
complex.) This can be done for example using the Metropolis method |17| . 
A random walk in the auxiliary-field space is simulated, such that the limit 



G 



distribution of sampled paths is proportional to |e lS \. This makes it possible 
to write (O) as a ratio of Monte Carlo averages, 

where {-) a .a' stands for the simulation average over auxiliary-field paths. A test 
of this approach has been carried out, showing that the term e ~ lRcS gener- 
ally makes the (averaged) numerator and denominator of Eq. (^) exceedingly 
small. Unless this "sign" problem can be overcome, the standard Metropolis 
method seems therefore to be inefficient in this context Since the weight of 
the paths in Eq. (|ll|) is complex (this is at the heart of the sign problem) a 
more promising possibility is the recourse to a simulation based on the complex 
Langevin equation Jll| |l9). In the Langevin algorithm (see, e.g., |2^]), paths 
distributed according to the "complex probability distribution" ~ e~ lS can be 
generated, allowing the computation of Eq. (|l3|) as a time average over a guided 
random walk for the fields in complex plane. In the case of interest here, the 
random walk for field a g is the solution of the stochastic differential equation 

da a i dS , . ,„ 
^ = -2^ + ^ (17) 

where t is a fictitious time (simulation time) and rj g is a (real) Gaussian white 
noise satisfying (r) g (t)) = and {ri g {t)ri g (t')) = — The first term in 
the right-hand side of Eq. ( p"7j ) can be seen as a "string" force which keeps a g 
close to the value for which the action iS is minimum, while the "noise" term 
is responsible for the sampling of a region in auxiliary-field space around this 
extremum. Although a general proof of the convergence of the complex Langevin 
simulation does not exist po|| , it turns out to work very nicely for a number of 
systems (the convergence is related to the location of the repulsive points of the 
Langevin dynamics). The Langevin simulation yields then a stochastic estimate 
of the output of the quantum computer, 

(0)~±J* +T dtO[a(t),*'(t)} (18) 

which is calculated by averaging 0[a, a'] for a sufficiently long random walk. 
Using Eq. (O) , the time derivative of the field a g can be written explicitly as 



da, 
~dt 



+ Ma,a'}+ Vg (t) (19) 



3 This numerical test has been performed on a small quantum circuit (L = 3, G = 4) using a 
one-field per gate Hubbard-Stratonovich transformation, but we are confident that the "sign" 
problem is generic. 
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with 

p \, /lA M^Mg^M^JQi) ram 



2=1 



One single term (I = a g ) differs from one in this product as only the 1-bit gate 
acting on qubit a g depends on o~ g [see Eqs. (0) and (|)]. One has 

= 11^9' (<V>V) ( 21 ) 



a 



with 



f T Mf„ r \- I -i a 9 A 9 U 9 as] (<r g >T- g ) if cf = g , > 

u g' Wg'^g'J-S TT [a„], s . l 22 J 

I Ug, {<Tgi,Tgi) otherwise 

The calculation of the derivative da g /dt (necessary to increment the fields along 
the random walk) requires thus the estimate of R g [a, a'] which is of the same 
kind as expression ( |l5| ) for the observable O: rather than inserting the observable 
O, one inserts the operator A g (conjugate to the field o~ g ) at a specific point 
in the ordered product of propagators. The coupling in the time evolution of 
the fields is obvious from Eq. (jTs|) . In particular, each pair of fields (a g ,T q ) 
is strongly coupled through the first term in the right-hand side of Eq. (|l9|). 
Indeed, combining Eq. (|l^) and its counterpart for r 9 , it is easy to see that the 
time-evolution of the field a g (or, equivalently, r g ) is governed by a second-order 
differential equation of the type d 2 a g /dt 2 ~ —a 2 a g /A, supplemented with a drift 
term (R g ) and a noise term (r) g ) in both the field a g and its velocity da g /dt. 

The detail of the Monte Carlo algorithm for implementing the complex 
Langevin simulation will be reported elsewhere. In short, the Langevin al- 
gorithm proceeds essentially in two alternating steps: (i) for the current value 
of the fields, calculate O and store it; (ii) update the fields by calculating all 
their time-derivatives da g /dt, using expression (|2^) for estimating the i? g 's. 
The time-average of O then yields the output of the quantum computation, the 
statistics being controlled by adjusting the length of the random walk. 



5 Example: the quantum FFT circuit 

Let us now discuss the scaling of the computational effort required to simulate 
a quantum algorithm, focusing on the quantum FFT algorithm jl4| used in 
Shor's factoring algorithm as an illustration.^] Since Shor's algorithm has been 
described in details in the literature (see, e.g., § for a review), it will be suffi- 
cient to note that, after a certain number of computational steps, the quantum 
register is in a periodic superposition of states \a) labeled by an integer between 

4 Note that the FFT is not the most computationally demanding task in Shor's algorithm, 
but this is unimportant for our illustrative purpose here. 
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Figure 2: Two-bit quantum Fast Fourier Transform circuit. It requires two 1-bit 
Hadamard gates and one 2-bit controlled-phase gate. 
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and 2 —1, the period being related to the sought factor of the composite 
number. The register is then subjected to a quantum FFT, resulting in a proba- 
bilistic estimate of the period (the probability of success can be made arbitrarily 
close to one by repeating the computation). The time-demanding task in the 
Monte Carlo simulation of the quantum FFT is the update of the AG auxiliary 
fields. Performing one step of the random walk in auxiliary-field space needs the 
computation of AG time-derivatives, requiring each the calculation of a single 
[cf. Eqs. ( pl| ) and (p2|)] one-bit matrix element (involving a product of about 
G/L non-unit 2x2 matrices). Thus, since G scales as L 2 /2 for the quantum 
FFT circuit, of the order of L 3 computation steps (2x2 matrix multiplication) 
are necessary to perform one step of the random walk. Assuming that the num- 
ber of steps necessary to achieve a given statistical error in the estimate of (O) 
does not grow exponentially with G (the sign problem should be overcome and 
the auto-correlation time of the random walk should not be exponential in G), 
the total number of computation steps would be polynomial in L. This does not 
rule out the possibility that, for a general quantum algorithm, the simulation 
effort might be polynomial in L whenever the number of gates G required in the 
quantum circuit is polynomial in L. This is an open question. 

As an example, we consider here a two-bit quantum FFT, i.e., the quantum 
computation of the discrete Fourier transform of a 4-point function (see Fig. ||) . 
The input qubits of the quantum register (L = 2) are labeled (and 1) for 
the least (and most) significant qubit. The two-bit quantum FFT circuit juj 
requires a single two-bit gate, a controlled-phase operator Coi — e lulAB acting 
on qubits a and b, with lu = n/2 and A = B = (1 — er z )/2, and two additional 
one-bit gates Hq and Hi, with H being the Hadamard transformation, 



The total unitary transformation is the ordered product U = HqCqxH\. The two 
one-bit gates Hq and H\ can be incorporated into the two-bit gate, which can 
in turn be written in terms of field-dependent one-bit gates using the Hubbard- 



JT|0) 
H\l) 



(|0) + |1»/V2 




(23) 
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Stratonovich representation, yielding U[a] — Uo(a)Ui(r), with 

Ui(r) = i=( e i r _ e L) (24) 

For a simple test of the Langevin algorithm, we consider here the Fourier 
transform of a constant function, i.e. the initial state is the product state 
2~ 1/2 (|0) + |1)) <g> 2~ 1 / 2 (|0) + |1)). The complex action can then be simply 
expressed as 

(1 4- p Mct-<0 \ 
-^—-2 ' (25) 

depending on the 4 auxiliary fields a, r, a', and r'. A straightforward calculation 
shows that the stochastic differential equations obeyed by the fields are 



= i 



i -( T -l/2)--tan(-(<7-<7')J +Va 
»2^-V2)+ 4tan(-(<7-<7 , )J+^ , 



da 
dc/_ 

IF 

dr u 

^ = (26) 

The Monte Carlo simulation of these equations is easy to perform. (Note that 
the fixed point of the Langevin dynamics — i.e., the path of minimum action — 
is a = a' = 0, r = r' = 1/2.) The resulting Monte Carlo averages for the 
1-bit observables Oo = |0)(0| and 0\ = |0)(0| converge to 1, implying that the 
expectation value for the output register is |00), as expected (the spectrum has 
a continuous component only). The simulation of larger quantum circuits using 
this technique is the subject of further work to be reported elsewhere. 



6 Conclusion 

We have shown that a quantum computer can be treated as a genuine quan- 
tum many-particle system, and that this approach sheds new light on quantum 
computation. More specifically, the use of a quantum Monte Carlo method 
might be interesting when considering "large" quantum computers because of 
the polynomial scaling of the auxiliary-field space in the dimension of the prob- 
lem. This advantage, however, is conditional on a circumvention of the Monte 
Carlo "sign" problem. In this respect, the use of a Langevin algorithm as a 
possibly efficient simulation technique is discussed. The stochastic simulation 
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of quantum computation proposed here could be useful for at least two reasons: 
(i) it could help in devising actual quantum computers by avoiding the need for 
an explicit experimental realization to test a quantum algorithm; (ii) it could 
give rise to a new class of "quantum-inspired" algorithms that could be imple- 
mented on an ordinary classical computer for solving certain computationally 
hard problems. 

We acknowledge C. Adami for many helpful discussions. This work has been 
funded in part by the NSF under Grant Nos. PHY 94-12818 and PHY 94-20470, 
and by a grant from DARPA/ARO through the QUIC Program (#DAAH04- 
96-1-3086). 
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